Method of Coalescence Microseismic Mapping Including Model&#39;s Uncertainty

ABSTRACT

A technique facilitates improved data acquisition and analysis with downhole tools and systems. The downhole tools and systems utilize arrays of sensing devices which are deployed in arrangements for improved sensing of data related to environmental features and/or tool parameters of tools located downhole in a borehole. For example, the tools and sensing systems may be operated to effectively sense and store characteristics related to components of downhole tools as well as formation parameters at, for example, elevated temperatures and pressures. Similarly, chemicals and chemical properties of interest in oilfield exploration and development also may be detected, measured and stored for analysis.

CROSS-REFERENCE TO RELATED APPLICATIONS

The present document is based on and claims priority to U.S. ProvisionalApplication Ser. No. 62/033,075, Method of Coalescence MicroseismicMapping Including Model's Uncertainty, filed Aug. 4, 2014, which isincorporated herein by reference in its entirety.

BACKGROUND

Hydrocarbon fluids such as oil and natural gas are obtained from asubterranean geologic formation, referred to as a reservoir, by drillinga well that penetrates the hydrocarbon-bearing formation. Analysis ofthe surrounding formation to find areas of interest (e.g., such ashydrocarbon containing deposits) and to avoid potential hazards are doneby a variety of methods. One type of analysis method involves the use ofacoustic signals (such as seismic, sonic, and microseismic) generated bya source (either actively or passively) and received by receivers.Complex mathematical manipulation can use factors such as the traveltimes of received waveforms emitted from the source to locate anddetermine various parameters and compositions of the geologicalcross-section.

SUMMARY

In general, a methodology and system are provided to facilitate improveddata acquisition and analysis with downhole tools and systems. Thedownhole tools and systems utilize arrays of sensing devices which aredeployed in arrangements for sensing of acoustic data related toenvironmental features and/or tool parameters of tools located downholein a borehole. For example, the tools and sensing systems may beoperated to effectively sense and store characteristics related tocomponents of downhole tools as well as formation parameters at, forexample, elevated temperatures and pressures. Similarly, chemicals andchemical properties of interest in oilfield exploration and developmentalso may be detected, measured and stored for analysis. The presentmethodology further enables improved mapping and analysis of theacoustic data.

However, many modifications are possible without materially departingfrom the teachings of this disclosure. Accordingly, such modificationsare intended to be included within the scope of this disclosure asdefined in the claims.

BRIEF DESCRIPTION OF THE DRAWINGS

Certain embodiments of the disclosure will hereafter be described withreference to the accompanying drawings, wherein like reference numeralsdenote like elements. It should be understood, however, that theaccompanying figures illustrate the various implementations describedherein and are not meant to limit the scope of various technologiesdescribed herein, and:

FIG. 1 is a graphical illustration of acceleration versus time withrespect to a noise window and a signal window, according to anembodiment of the disclosure;

FIG. 2 is a graphical illustration showing waveforms containing P- andS-event signals, according to an embodiment of the disclosure;

FIG. 3 is a graphical illustration showing a characteristic function,according to an embodiment of the disclosure;

FIG. 4 is a graphical illustration showing an unmodified characteristicfunction, according to an embodiment of the disclosure;

FIG. 5 is a graphical illustration showing a modified characteristicfunction, according to an embodiment of the disclosure;

FIG. 6 is a graphical illustration showing an example of a plurality ofsynthetic waveforms, e.g. synthetic waveforms to test an onset timeestimation, according to an embodiment of the disclosure;

FIG. 7, is a graphical illustration showing location error versusevents, according to an embodiment of the disclosure;

FIG. 8, is a schematic illustration of a downhole system utilizing a raytracing for a given location and a downhole array, according to anembodiment of the disclosure; and

FIG. 9, is a graphical illustration showing RMS errors of onset timeestimation versus events, according to an embodiment of the disclosure.

DETAILED DESCRIPTION

Reference throughout the specification to “one embodiment,” “anembodiment,” “some embodiments,” “one aspect,” “an aspect,” or “someaspects” generally indicates a particular feature, structure, method, orcharacteristic described in connection with the embodiment or aspect isincluded in at least one embodiment of the present disclosure. Thus, theappearance of the phrases “in one embodiment” or “in an embodiment” or“in some embodiments” in various places throughout the specification arenot necessarily referring to the same embodiment. Furthermore, theparticular features, structures, methods, or characteristics may becombined in a suitable manner in one or more embodiments. The words“including” and “having” generally have the same meaning as the word“comprising.”

As used throughout the specification and claims, the term “downhole”refers to a subterranean environment, e.g. a subterranean environmentalong a wellbore. “Downhole tool” is used broadly to refer to a toolused in a subterranean environment including, but not limited to, alogging tool, an imaging tool, an acoustic tool, a permanent monitoringtool, and a combination tool. Additionally, “coalesce” or “coalescence”may be used to refer to a coming together of different measures into amap such as a spatial map. The terms also may be used to indicateevaluating the time evolving contributions from multiple sensors orcomponents without including non-linear event location methods.Furthermore, the term “sensor” or “receiver” may be used to refer to asingle device location which may have one or more “detectors” able toreceive measurements at that location.

The various techniques disclosed herein may be utilized to facilitateand improve data acquisition and analysis in downhole tools and systems.Examples of downhole tools and systems utilize arrays of sensingdevices, e.g. receivers, configured or designed for easy attachment anddetachment in downhole sensor tools or modules that are deployed forpurposes of sensing data relating to environmental and tool parametersdownhole, e.g. within a borehole. The sensing systems disclosed hereinmay effectively sense and store characteristics relating to componentsof downhole tools as well as formation parameters, e.g. formationparameters involving elevated temperatures and pressures. Chemicals,chemical properties, and a variety of other subterranean features ofinterest in oilfield exploration and development also may be measured,stored, and evaluated via data obtained from the sensing systems andprocessed according to data processing techniques described herein.

The sensing system embodiments described herein may be incorporated intool systems such as wireline logging tools, measurement-while-drillingtools, logging-while-drilling tools, permanent monitoring systems, drillbits, drill collars, sondes, among others. For purposes of thedescription herein, when the terms wireline, cable line, slickline,coiled tubing, or conveyance are used it should be understood that theseand other suitable conveyance systems and techniques may be employed. Asdescribed in greater detail below, certain embodiments provide a methodof coalescence microseismic mapping which include, e.g. account for, amodel's uncertainty. In some cases, the results of the methods anddescriptions contained herein may be used for determining hypocenterlocations with single well observation during hydraulic fracturingmonitoring, indicating the growth and propagation of fractures. Thisinformation can be used for better control of fracturing operation andto ensure that critical water deposits are not breached duringfracturing operation.

To apply an inverse method using probability density functions (see, forexample, Tarantola and Valette: Inverse Problems=Quest for Information.J. Geophys. 42, 159-170, 1982), maximum likelihood, and CoalescenceMicroseismic Mapping (CMM) (see, for example, Drew et al., Automatedmicroseismic event detection and location by continuous spatial mapping:Proceedings of the Annual Technical Conference and Exhibition, SPE,95113, 2005) methods to locate hypocenters with insufficient surveyconfigurations, the polarizations of the P- and S-waves may be used tosolve the directional ambiguities. Because the qualities of estimatedpolarizations are substantially varied over receivers, a simple averageof directions estimated using polarizations causes substantial errors ofthe locations.

The error of polarization is defined by the difference between themeasured and computed polarizations for the possible hypocenterlocations. The cost function is defined using the errors ofpolarizations and the weight functions computed using signal qualities(e.g., SNR, linearity and orthogonality). The minimization problem maybe solved by a grid search in the 3-D space to obtain the most probablelocation or direction of the hypocenter. From the estimated direction,the probability density function of polarization may be formulated andcombined with the marginal probability density function obtained usingthe picked event times (or product of SNR in CMM) to locate thehypocenters.

To locate hypocenters in a 3-D space using the observed data andvelocity model, the information from polarization is used whensufficient coverage of the hypocenter cannot be obtained. This is thecase for a single well measurement. Because the qualities of theestimated polarization are substantially varied over receivers, a simpleaverage of directions estimated using polarizations causes substantialerrors in determining the locations. Although attempts have been made toestablish a relation between linearity and uncertainty of polarizationmeasurements by using synthetic data (see, for example, Drew et al.,Microseismic event azimuth estimation: establishing a relationshipbetween hodogram linearity and uncertainty in event azimuth; 2008 SEGTechnical Program Expanded Abstracts, 1446-1449), the results are notapplicable to real datasets due to the fact that the noises in real dataare coherent and their linearity is very high. Embodiments describedherein provide a new method to estimate the direction of hypocenter frommeasured polarizations.

According to an embodiment, a hypocenter location method using simplytravel-times is run first to obtain a start time and hypocenter. Thetables of travel times and polarizations at the receivers (e.g. look-uptables) may be prepared for the possible hypocenter locations as grids.The travel times and polarizations may be computed using ray tracing, anEikonal equation or a finite difference method. The polarization at areceiver associated with a picked travel time may be compared to thepolarizations stored in the look-up table. The cost function may then becomputed using the difference between the measured and computedpolarizations for the subject receivers. By minimizing the cost functionthrough the use of the grid search, the most probable location ordirection of the hypocenter can be determined. The processed dataacquired by minimizing the cost function to determine probable locationsand directions of hypocenters may be used via, for example, aprocessor-based system to evaluate subterranean features, e.g. featuresrelated to downhole tools, fractures, geologic or other environmentalparameters, petroleum reservoir features, and/or other desired features.

In a first stage of an embodiment, an inverse method using a probabilitydensity function (e.g. see Tarantola and Valette, 1982) (or CMM) isapplied to find potential candidates for hypocenter location and thestart time of the event, then, onset times of event are estimated fromthe potential candidates. Because a hypocenter location is computedusing simply the travel times with insufficient survey configurations,there are directional ambiguities for hypocenter locations where thestart time of an event cannot be uniquely determined Regarding thesignal components to be used in the computation of direction ofhypocenter location, both the P- and S-waves, the P-waves alone, or theS-waves alone may be specified by the users. The obtained start time isused to restrict grid points in the grid search. The computed event timefor the i-th receiver may be defined by:

T _(i)(x,y,z)=T ₀ +t _(i)(x,y,z)

where T_(i)(x,y,z) is a computed onset time of an event signal at thei-th receiver for the hypocenter at (x,y,z) and T₀ is the start time ofevent. t_(i)(x,y,z) is the travel time between the hypocenter (x,y,z)and the i-th receiver.

The grid search in the following procedure is restricted for the gridpoints which satisfy:

|T _(i)(x,y,z)−Pick_(i) |<ΔT

where Pick_(i) is the picked time of event signal and ΔT is a parameterdefined by the users.

A look-up table, which stores the travel-times from hypocenters to areceiver and the polarization vectors at a receiver, may be created in avolume where hypocenters are anticipated. By way of example, a look-uptable may be created using ray tracing, an Eikonal equation, or a finitedifference method. The error of polarization vectors may be defined by:

θ_(i,Error)(x,y,z)=|A cos(P _(i,mes) ·P _(i,comp)(x,y,z))|

where i is the index of a receiver, θ_(i,Error) is the polarizationerror, and P_(i,mes) and P_(i,comp)(x,y,z) are the measured and computedpolarization vectors at the i-th receiver for the hypocenter location at(x,y,z), respectively. The vectors used in this example are normalized.

The hypocenter location or most probable direction of the hypocenter maybe solved by minimizing the cost function defined by:

${\theta_{Error}\left( {x,y,z} \right)} = \left( \frac{\sum\limits_{i = 1}^{N}{w_{i}{\theta_{i,{Error}}^{n}\left( {x,y,z} \right)}}}{\sum\limits_{i = 1}^{N}w_{i}} \right)^{- n}$

where n is an appropriate power to be selected and w_(i) is a weightfunction of signal-to-noise (SNR), linearity, and orthogonality of anevent signal denoted by:

w _(i)=ƒ(SNR_(i),Linearity_(i),Orthogonality_(i)).

The weight function may be simply written as:

w _(i)=SNR_(i) ^(p)·Linearity_(i) ^(q)·Orthogonality_(i) ^(r)

where p, q and r are appropriate powers to be selected. SNR in the aboveequation may be replaced by:

SNR′=max(SNR−1,0).

By denoting (x′,y′,z′) that minimizes θ_(Error)(x,y,z), the jointprobability density function (PDF) of polarization to be applied in themethod using probability density function may be defined by:

${{PDF}\left( {x,y,z} \right)} = {\prod\limits_{i}{\exp \left\lbrack {{- \frac{1}{2}}\left( \frac{{\theta_{i,{Error}}\left( {x,y,z} \right)} - {\theta_{i,{Error}}\left( {x^{\prime},y^{\prime},z^{\prime}} \right)}}{\sigma_{i}} \right)^{2}} \right\rbrack}}$where σ_(i) = max (θ_(i, Error)(x^(′), y^(′), z^(′)), θ_(min))

and θ_(min) is the minimum angle defined by the users.

By defining a standard deviation as:

$\sigma_{\theta} = \left( {\sum\limits_{i}\sigma_{i}^{- 2}} \right)^{- 2}$

in the vicinity of (x′,y′,z′)

$\frac{{\theta_{Error}\left( {x,y,z} \right)} - {\theta_{Error}\left( {x^{\prime},y^{\prime},z^{\prime}} \right)}}{\sigma_{\theta}} \approx {\sum\limits_{i}\left( \frac{\begin{matrix}{{\theta_{i,{Error}}\left( {x,y,z} \right)} -} \\{\theta_{i,{Error}}\left( {x^{\prime},y^{\prime},z^{\prime}} \right)}\end{matrix}}{\sigma_{\theta}} \right)^{2}}$

is anticipated for rather good quality data. Therefore the joint PDF maybe defined by:

${{PDF}\left( {x,y,z} \right)} = {\exp \left\lbrack {{- \frac{1}{2}}\left( \frac{{\theta_{Error}\left( {x,y,z} \right)} - {\theta_{Error}\left( {x^{\prime},y^{\prime},z^{\prime}} \right)}}{\sigma_{\theta}} \right)^{2}} \right\rbrack}$

or more simply it may be defined by:

${{PDF}\left( {x,y,z} \right)} = {\exp \left\lbrack {{- \frac{1}{2}}\left( \frac{\theta_{Error}\left( {x,y,z} \right)}{\sigma_{\theta}} \right)^{2}} \right\rbrack}$

By choosing appropriate n, a robust estimation for bad quality data isstably made using the above equation.

In some embodiments, the Coalescence Microseismic mapping method using amodel's uncertainty may be given by the following equation:

${F\left( {x,y,z} \right)} = {\int{\prod\limits_{i = 1}^{N}{{f_{i}\left( \tau_{i} \right)} {\exp\left( \begin{matrix}{{- \frac{1}{2}}\left( {\overset{\sim}{\tau} - {\overset{\sim}{h}\left( {x,y,z} \right)}} \right)^{T}} \\{C_{T}^{- 1}\left( {\overset{\sim}{\tau} - {\overset{\sim}{h}\left( {x,y,z} \right)}} \right)}\end{matrix} \right)}{\tau_{1}}\mspace{14mu} \ldots \mspace{14mu} {\tau_{N}}}}}$

The multiple integral may be applied for the time duration where eventsignals are found. The polarization PDF also may be applied using thepreviously described method.

In this example, the characteristic function is defined as:

${f(t)} = \sqrt{\frac{{nw}{\int_{t}^{t + {sw}}{{A^{2}(\tau)}\ {\tau}}}}{{sw}{\int_{t - {nw}}^{t}{{A^{2}(\tau)}\ {\tau}}}}}$

where ƒ(t) is the characteristic function, nw and sw are window lengths,A(t) is the envelope of seismic traces. A(t) is computed for the 3component data or the projected data based on the polarization estimatedby signal or model. A smoothing operator, e.g. a Hanning filter, may beapplied to A(t). (See the graphical representations of accelerationversus time and SNR versus time in FIGS. 1-3).

The characteristic function may be modified as:

g(t)=max(ƒ(t)−1,ε)

where ε is appositive number sufficiently smaller than one. Examples ofthe SNR versus time and the modified SNR versus time are illustratedgraphically in FIGS. 4 and 5, respectively.

By regarding the modified characteristic function as PDF, CMM containingthe model error may be defined as:

F(x,y,z)=∫Π_(i=1) ^(N)ƒ_(i)(τ_(i))exp(−½({tilde over (τ)}−{tilde over(h)}(x,y,z))^(τ) C _(T) ⁻¹({tilde over (τ)}−{tilde over (h)}(x,y,z)))dτ₁ . . . dτ _(N)

where i and N are the index and number of modified characteristicfunctions. A single index may be used when the multi-components areused. τ_(i) is the recording time, and ƒ_(i)(τ_(i)) is the modifiedcharacteristic function. {tilde over (τ)} is the vector defined as:

${\overset{\sim}{\tau}}_{i} = {\tau_{i} - \overset{\_}{\tau}}$$\overset{\_}{\tau} = \frac{\sum\limits_{i = 1}^{N}{w_{i}\tau_{i}}}{\sum\limits_{i = 1}^{N}w_{i}}$

where w_(i) is the weight defined as below.

In this example, {tilde over (h)}(x,y,z) is the vector and may bedefined as:

${\overset{\sim}{h_{i,j}}\left( {x,y,z} \right)} = {{h_{i}\left( {x,y,z} \right)} - {\overset{\_}{h}\left( {x,y,z} \right)}}$${\overset{\_}{h}\left( {x,y,z} \right)} = \frac{\sum\limits_{i = 1}^{N}{w_{i}{h_{i}\left( {x,y,z} \right)}}}{\sum\limits_{i = 1}^{N}w_{i}}$

where h_(i)(x,y,z) is the travel time. The weight may be computed usingthe covariance matrix as

$w_{j} = {\sum\limits_{i = 1}^{N}p_{ij}}$ C_(T) = σ_(ij)C_(T)⁻¹ = p_(ij)

where σ_(ij) is computed as:

$\sigma_{i,j} = {\sigma_{T}^{2}{\exp\left( {{- \frac{1}{2}}\frac{D_{ij}^{2}}{\Delta^{2}}} \right)}}$

σ_(τ) is the model's uncertainty specified by the users, D_(ij) is thedistance between the receiver positions of the i and j traces, and Δ isthe correlation length.

The multiple integral with respect to τ_(i) may be applied for the timeduration where event signals are found.

Synthetic data may be used to validate the methodology. According to anexample, a plurality of receivers, e.g. 12 receivers, may be installedin a vertical well and a 2-D (r-z) problem is solved. The Monte-Carlomethod may be used for the multiple time integrations. Referringgenerally to FIG. 6, a graphical illustration is provided of the 12synthetic waveforms. The data may be processed to determine locationerrors for the seismic events, as illustrated in FIG. 7.

In another embodiment, an onset estimation method using CoalescenceMicroseismic Mapping may be utilized. This method uses 1) Modifiedcharacteristic function, 2) Onset time estimation using the location andorigin time that are estimated by CMM, 3) temporal buffer lists to mergeand select detection results, and 4) event location method using theonset estimation results.

In at least some embodiments, a velocity model may be used to computethe travel-times and polarizations of, for example, the P-, SH- andSV-waves. The waves are selected depending on the survey configurationand the data condition (e.g. P-wave alone, P- and SH-waves combined, orother suitable waves/wave combinations). The grids in the 3-D space areused to estimate the onset times.

A maximum grid size may be defined as:

${\delta \; x} = {\frac{\sqrt{2}}{2}{\min \left( {{V_{P,\min}{sw}_{P}},{V_{S,\min}{sw}_{S}}} \right)}\mspace{14mu} {for}{\mspace{11mu} \;}2\text{-}D}$${\delta \; x} = {\frac{2\sqrt{3}}{3}{\min \left( {{V_{P,\min}{sw}_{P}},{V_{S,\min}{sw}_{S}}} \right)}\mspace{14mu} {for}\mspace{14mu} 3\text{-}D}$

where V_(p,min) and V_(s,min) are minimum P- and S-velocities; andsw_(p) and sw_(s) are signal windows for P- and S-waves respectively.The signal window is defined in greater detail below.

The projection of 3 component waveforms to a scalar function iseffective for extracting a targeted signal component (e.g., P-, SH- orSV-waves) and for avoiding contamination by the other signal components.The inner product is taken with the 3 component waveforms and referencevector as:

g(t,x,y,z)=ref(x,y,z)·(g _(x)(t),g _(y)(t),g _(z)(t))

where g(t, x, y,z) is the project waveform, ref(x,y,z) is the referencevector, and (g_(x)(t),g_(y)(t),g_(z)(t)) is the 3 component waveforms.ref(x,y,z) is the polarization direction of signal component to beextracted. ref(x,y,z) may be taken from a-priori information such as thedirection from the receivers to an injection point or the polarizationscomputed using the model (in this case, the projected waveforms arecomputed for the location). The azimuthal component of ref(x,y,z) may beused alone because the inclination information tends to be lessaccurate.

A characteristic function may be defined as:

${f(t)} = \sqrt{\frac{{nw}{\int_{t}^{t\rightarrow{sw}}{{A^{2}(\tau)}{\tau}}}}{{sw}{\int_{t - {nw}}^{t}{{A^{2}(\tau)}{\tau}}}}}$

where ƒ(t) is the characteristic function, nw and sw are window lengths,and A(t) is the envelope of seismic traces. A(t) is computed for the 3component data or the projected data based on the polarization estimatedby signal or model. A smoothing operator, such as a Hanning filter, maybe applied to A(t). ƒ(t) is also called with respect to SNR as it is theratio of root mean square (RMS)s of the window graphically illustratedin FIG. 1. FIG. 1 illustrates the signal window and the preceding noisewindow. Often sw may be chosen as a half of dominant period of signaland nw is chosen as three to six times of Fay. It should be noted FIGS.2 and 3 illustrate examples of the signal and characteristic functions.Additionally, FIGS. 4 and 5, respectively, represent examples of theunmodified and modified characteristic functions. As referenced above,the characteristic function may be modified as:

g(t)=max(ƒ(t)−1,ε)

where ε is a positive number sufficiently smaller than one. ε can be e⁻²or e⁻⁴.

A time-space CMM function may be defined as:

${F\left( {t,x,y,z} \right)} = {\prod\limits_{i = 1}^{N}\; {\prod\limits_{j = 1}^{M}\; {f_{ij}\left( {t + {T_{ij}\left( {x,y,z} \right)}} \right)}}}$or${F\left( {t,x,y,z} \right)} = {\min\limits_{j}\left\lbrack {\prod\limits_{i = 1}^{N}\; {f_{ij}\left( {t + {T_{ij}\left( {x,y,z} \right)}} \right)}} \right\rbrack}$

where (t,x,y,z) is the origin time and location to be tested, and areindices of receiver and wave component (e.g., P-, SH- and SV-waves), Nand M are numbers of receivers and wave components, ƒ_(ij)(t) is themodified characteristic function, and T_(ij)(x,y,z) is the travel-timeof the j-th component wave from the receiver to the location.

Components to be used are selected according to the data (e.g. P-wavesalone, P- and SH-waves, or other suitable waves/wave combinations). Theorthogonality of polarization may be optionally added as:

${F\left( {t,x,y,z} \right)} = {\prod\limits_{i = 1}^{N}\; {{{Orth}_{i}\left( {t,x,y,z} \right)}{\prod\limits_{j = 1}^{M}\; {f_{ij}\left( {t + {T_{ij}\left( {x,y,z} \right)}} \right)}}}}$${{Orth}_{i}\left( {t,x,y,z} \right)} = {\prod\limits_{j < k}\; \left\lbrack {1 - \left\{ {{p_{i}\left( {t + {T_{ij}\left( {x,y,z} \right)}} \right)} \cdot {p_{i}\left( {t + {T_{ik}\left( {x,y,z} \right)}} \right)}} \right\}^{2}} \right\rbrack}$

where p_(i)(t) is the polarization vector.

In FIG. 8, an embodiment of a monitoring system 30 is provided toillustrate ray tracing for a given formation location 32 which serves asthe source/event 32 for a plurality of rays 34 which are received by adownhole array 36. The downhole array 36 may comprise a plurality oftool/receivers 38 disposed in, for example, a wellbore 40. When thenumber of receivers 38 is large, the multiplications may be changed tosummation by taking logarithms. As shown, this embodiment may be asingle well monitoring system using a wireline or other conveyeddownhole array 36.

According to an embodiment, the threshold of time-space CMM function maybe given by:

Threshold_(Detection)=Threshold_(SNR) ^((1−α)NM)ε^(αNM)

If a strict detection option is used, the threshold may be given by:

Threshold_(Detection)=Threshold_(SNR) ^((1−α)N)ε^(αN)

where Threshold_(SNR), is the reference SNR value and a is the failureratio.

An orthogonality constraint is optionally applied to the time-space CMMfunction as follows:

Orth > OrthLimit${Orth} = {\frac{1}{\sum\limits_{i = 1}^{N}w_{i}}{\sum\limits_{i = 1}^{N}{w_{i}{{Orth}_{i}\left( {t,x,y,z} \right)}}}}$$w_{i} = \sqrt{{f_{ij}\left( \tau_{j} \right)}{f_{ik}\left( \tau_{k} \right)}}$

The summation may be taken for the orthogonality whose SNR is greaterthan a threshold as w_(i)>Threshold_(SNR).

Onset times may be estimated at a given time, KΔt, for a location thatsatisfies:

${\max\limits_{x,y,z}{F\left( {{k\; \Delta \; t},x,y,z} \right)}} \geq {Threshold}_{Detection}$or F(k Δ t, x, y, z) ≥ Threshold_(Detection)

For the former and latter cases, the onset times may be estimated forsingle and multiple grid points respectively. Suppose a point (x₀,y₀,z₀)satisfies the above equations, the onset times may be estimated usingtravel times, which, in turn, may be computed using a velocity model as:

=kΔt+T _(ij)(x ₀ ,y ₀ ,z ₀)

where

is the estimated onset time for the j-th signal component of the i-threceiver, k is the index of time, Δt is the sampling period, andT_(ij)(x₀,y₀,z₀) is the travel time of j-th signal component between thelocation (x₀,y₀,z₀) and the i-th receiver position. For the signals thatsatisfy ƒ_(ij)(

)≧Threshold_(SNR,ij), the accurate onset times are determined by findingthe maximum of ƒ_(ij)(t) in a window centered by

. The correction value δτ(−sw_(ij)<δτ<sw_(ij)) is selected so that

${f_{ij}\left( {+ {\delta \; \tau}} \right)} = {\max\limits_{{- {sw}_{ij}} < {\delta \; t} < {sw}_{ij}}{f_{ij}\left( {+ {\delta \; t}} \right)}}$

The signal window length of the characteristic function may be used forthe length of window, sw_(ij). If a smoothing operator is applied toƒ_(ij)(t), the length of window sw_(ij) is the addition of the signalwindow length and operator length. Then the onset time may be obtainedas:

τ_(ij)=

+δτ

for the signal which satisfies:

ƒ_(ij)(τ_(ij))≧Threshold_(SNR,ij)

where Threshold_(SNR,ij) is the threshold for SNRs.

For the obtained onset times, the orthogonality constraint may beapplied for the multi component case. If the constraint is notsatisfied, the onset times are rejected. Using the characteristicfunction for the obtained onset times, a new detection function may becomputed as:

${Detect} = {\prod\limits_{ij}\; {\overset{\_}{f_{ij}}\left( \tau_{ij} \right)}}$or${Detect} = {\min\limits_{j}\left\lbrack {\prod\limits_{i}\; {\overset{\_}{f_{ij}}\left( \tau_{ij} \right)}} \right\rbrack}$where${\overset{\_}{f_{ij}}\left( \tau_{ij} \right)} = {{{f_{ij}\left( \tau_{ij} \right)}\mspace{14mu} {if}\mspace{14mu} {f_{ij}\left( \tau_{ij} \right)}} \geq {{Threshold}_{{SNR},{ij}}\mspace{14mu} \left( {{for}\mspace{14mu} {all}\mspace{14mu} j} \right)}}$${\overset{\_}{f_{ij}}\left( \tau_{ij} \right)} = {{ɛ\mspace{14mu} {if}\mspace{11mu} {f_{ij}\left( \tau_{ij} \right)}} < {{Threshold}_{{SNR},{ij}}\mspace{14mu} \left( {{for}\mspace{14mu} {all}\mspace{14mu} j} \right)}}$

If the number of onset times that satisfy:

ƒ_(ij)(τ_(ij))≧Threshold_(SNR,ij)(for all j)

is less than the threshold, Threshold_(N), the selected onset times arecanceled.

Lists may be established to save onset times and other attributes intemporal lists as follows:

Onset_(ij)[l] is the l-th entry of τ_(ij), SNR_(ij)[l] is ƒ_(ij)(τ_(ij)), Detection[l] is Detect, and M is the number of entries. WhenDetect≧Threshold_(Detection), the new entries are saved.

Various options may be used to operate temporal lists. According to afirst embodiment, an option is to replace overlapped events and compresstemporal lists. If, for example, Onset_(ij)[l] and Onset_(ij)[m] (l<m)for any i and all j satisfy:

SNR_(ij) [l]≧Threshold_(SNR,ij)

SNR_(ij) [m]≧Threshold_(SNR,ij)

|Onset_(ij) [l]−Onset_(ij) [m]≦sw _(ij)

and

Detection[m]>Detection[l]

Then the m-th entries are copied to the l-th entries and the m-thentries are removed from the lists.

According to a second embodiment, another option is to merge overlappedevents and compress temporal lists. If, for example, Onset_(ij)[l] andOnset_(ij)[m] for any i and all j satisfy:

SNR_(ij) [l]≧Threshold_(SNR,ij)

SNR_(ij) [m]≧Threshold_(SNR,ij)

|Onset_(ij) [l]−Onset_(ij) [m]≦sw _(ij)

Then for n that satisfies:

SNR_(ij) [l]<Threshold_(SNR,ij)

SNR_(ij) [m]≧Threshold_(SNR,ij)

entries of the list are copied as:

Onset_(ij) [l]=Onset_(ij) [m]

SNR_(ij) [l]=SNR_(ij) [m]

and Detection[l] is computed for new SNR_(ij)[l]. Ultimately, the m-thentries may be removed from the lists. This operation may be continueduntil the various l and m combinations are tested. If inconsistencyhappens, the event may be saved as a new event.

Another option may be to merge distinct events in a certain time windowand compress temporal lists. If, for example, any i and n (i≠n) for allj satisfy:

SNR_(ij)(l) ≥ Threshold_(SNR, ij)  and  SNR_(ij)[m] < Threshold_(SNR, ij)SNR_(nj)(l) < Threshold_(SNR, nj)  and  SNR_(nj)[m] ≥ Threshold_(SNR, nj)${\max\limits_{i,n,j}\left( {{{Onset}_{ij}\lbrack l\rbrack} - {{Onset}_{nj}\lbrack m\rbrack}} \right)} \leq T$

entries of lists are copied as:

Onset_(nj) [l]=Onset_(ij) [m]

SNR_(nj) [l]=SNR_(ij) [m]

and Detection[l] is computed for new SNR_(ij)[l] where T is a thresholdspecified by the users.

Ultimately, the m-th entries are removed from the lists. It should benoted that indices i which satisfy

SNR_(ij) [l]<Threshold_(SNR,ij) and SNR_(ij) [m]<Threshold_(SNR,ij)

can be ignored. This operation may be continued until all l and mcombinations are tested. If inconsistency happens, the event may besaved as a new event.

Some embodiments involve removal of multiple detections for a singleonset time. There are cases in which an onset time is detected as the P-and SH-waves of different events. In such a case, the orthogonality maybe used to select events. In this example, one event is saved in thelists and others are removed.

Onset times also may be provided using the time-space CMM function.Accordingly, the event location method by Tarantola & Valette (1982) maybe applied in some embodiments. In this example, a single index may beused for the receivers and signal components according to the following:

${\sigma_{p}\left( {x,y,z} \right)} = {\exp \left( {{- \frac{1}{2}}\left( {\overset{\sim}{\tau} - {\overset{\sim}{h}\left( {x,y,z} \right)}} \right)^{T}{P\left( {\overset{\sim}{\tau} - {\overset{\sim}{h}\left( {x,y,z} \right)}} \right)}} \right)}$where P = (C_(t) + C_(T))⁻¹ = [p_(ij)]${\left( {x,y,z} \right)} = {{h_{i}\left( {x,y,z} \right)} - {\overset{\_}{h}\left( {x,y,z} \right)}}$${\overset{\_}{h}\left( {x,y,z} \right)} = \frac{\sum\limits_{i = 1}^{N}{w_{i}{h_{i}\left( {x,y,z} \right)}}}{\sum\limits_{i = 1}^{N}w_{i}}$$w_{j} = {\sum\limits_{i = 1}^{N}p_{ij}}$

The covariance matrix C_(t) is diagonal and its elements may be:

for  f_(i)(τ_(i)) > γ ⁻²$\sigma_{ii} = \frac{sw}{\sqrt{{- 2}\; \log \frac{\gamma \; ^{- 2}}{f_{i}\left( \tau_{i} \right)}}}$else σ_(ii) = sw$\sigma_{ij} = {\sigma_{T}^{2}{\exp \left( {{- \frac{1}{2}}\frac{D_{ij}^{2}}{\Delta^{2}}} \right)}}$

As with certain embodiments described above, a synthetic test may beconducted to validate the methodology. For example, synthetic waveformsmay be generated for a plurality of level downhole tools, e.g. 12downhole tools/receivers, using a constant velocity model. Random noisesmay be added. For the onset time estimation, the P- and S-velocities maybe perturbed by a desired amount, e.g. +2% and −2%, respectively.Synthetic waveforms, such as synthetic waveforms similar to thoseillustrated in FIG. 6 enable testing of the onset time estimation. InFIG. 9, a graphical example is provided illustrating the RMS errorsversus events. The graphical example shows RMS errors of onset timeestimation. In this example, the curves or lines 42, 44 show the RMSerrors of the P- and SH-waves respectively. It should be noted theevaluation of acoustic data received by receivers 38, the establishmentof grids, the determination/use of cost functions, the use offormulas/algorithms described herein, and/or the output of informationregarding subterranean features may be carried out automatically via aprocessor-based system, e.g. a processor-based computer.

Although a few embodiments of the disclosure have been described indetail above, those of ordinary skill in the art will readily appreciatethat many modifications are possible without materially departing fromthe teachings of this disclosure. Accordingly, such modifications areintended to be included within the scope of this disclosure as definedin the claims.

What is claimed is:
 1. A method for evaluating subterranean features,comprising: obtaining onset times of event using the candidates of starttime and hypocenter given by Coalescence Microseismic Mapping (CMM);computing travel times and polarizations of acoustic waves using raytracing; preparing tables of travel times and polarizations at aplurality of receivers based on the ray tracing to establish hypocenterlocations as a grid; comparing the polarizations associated with traveltimes of the acoustic waves received at receivers with polarizationspreviously computed and stored; determining a cost function using thedifferences between measured at obtained onset times and computedpolarizations at the receivers; and establishing probable locations ofhypocenters by minimizing the cost function via grid searching.
 2. Themethod as recited in claim 1, wherein establishing further comprisesestablishing probable directions of the hypocenters.
 3. The method asrecited in claim 1, wherein using comprises applying an inverse methodusing a probability function to find potential candidates for hypocenterlocations and start times of events.
 4. The method as recited in claim1, wherein using comprises using acoustic waves in the form of P-waves.5. The method as recited in claim 1, wherein using comprises usingacoustic waves in the form of S-waves.
 6. The method as recited in claim1, wherein using comprises using acoustic waves in the form of P-wavesand S-waves.
 7. The method as recited in claim 1, wherein determiningthe cost function comprises determining the cost function defined by:${\theta_{Error}\left( {x,y,z} \right)} = \left( \frac{\sum\limits_{i = 1}^{N}{w_{i}{\theta_{i,{Error}}^{n}\left( {x,y,z} \right)}}}{\sum\limits_{i = 1}^{N}w_{i}} \right)^{- n}$where n is a power selected and w_(i) is a weight function ofsignal-to-noise (SNR), linearity, and orthogonality of an event signaldenoted by:w _(i)=ƒ(SNR_(i),Linearity_(i),Orthogonality_(i)).
 8. The method asrecited in claim 1, wherein computing travel times and polarizationscomprises using a velocity model.
 9. The method as recited in claim 1,further comprising estimating onset times for single and multiple gridpoints of the grid.
 10. The method as recited in claim 9, furthercomprising saving the onset times in temporal lists.
 11. A method,comprising: deploying an array of sensing devices to monitor acousticwaves resulting from acoustic events; measuring travel times andpolarizations of the acoustic waves received by the array of sensingdevices; computing a grid of stored travel times and polarizations;establishing a cost function using differences between measured andcomputed travel times and polarization; minimizing the cost function viaa grid search to establish probable locations and directions ofhypocenters; and outputting information related to the hypocenters. 12.The method as recited in claim 11, wherein deploying comprises deployingthe array of sensing devices in a wellbore.
 13. The method as recited inclaim 11, wherein measuring comprises measuring P-waves.
 14. The methodas recited in claim 11, wherein measuring comprises measuring S-waves.15. The method as recited in claim 11, wherein establishing comprisesestablishing the cost function by:${\theta_{Error}\left( {x,y,z} \right)} = \left( \frac{\sum\limits_{i = 1}^{N}{w_{i}{\theta_{i,{Error}}^{n}\left( {x,y,z} \right)}}}{\sum\limits_{i = 1}^{N}w_{i}} \right)^{- n}$where n is a power selected and w_(i) is a weight function ofsignal-to-noise (SNR), linearity, and orthogonality of an event signaldenoted by:w _(i)=ƒ(SNR_(i),Linearity_(i),Orthogonality_(i)).
 16. The method asrecited in claim 11, wherein computing comprises computing via aprocessor-based computer.
 17. The method as recited in claim 11, furthercomprising using data obtained by minimizing the cost function toevaluate a subterranean feature.
 18. A method, comprising: comparing ameasured polarization, associated with a travel time of an acousticsignal received at a receiver, with a stored polarization; computing adifference between the measured polarization and the storedpolarization; using the difference to determine a probable location of ahypocenter; and outputting data on the probable location to facilitateevaluation of a subterranean feature.
 19. The method as recited in claim18, wherein using comprises determining a cost function to evaluate thedifference between the measured polarization and the storedpolarization.
 20. The method as recited in claim 18, wherein comparingcomprises comparing measured polarizations received at a plurality ofreceivers arranged in a subterranean array.